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■ Abstract 

A new method to determine the low-energy couphngs of the AS" = 1 weak Hamihonian is 
presented. It rehes on a matching of the topological poles in 1/m^ of three-point correlators of 
two pseudoscalar densities and a four-fermion operator, measured in lattice QCD, to the same 
observables computed in the e-regime of chiral perturbation theory. We test this method in a 
theory with a light charm quark, i.e. with an SU(4) flavour symmetry. Quenched numerical 
measurements are performed in a 2 fm box, and chiral perturbation theory predictions are 
worked out up to next-to-leading order. The matching of the two sides allows to determine 
the weak low-energy couplings in the SU(4) limit. We compare the results with a previous 
determination, based on three-point correlators containing two left-handed currents, and 
discuss the merits and drawbacks of the two procedures. 



1. Introduction 



Understanding quantitatively or even just qualitatively non-Ieptonic kaon decays, K vrvr, 
starting from the QCD Lagrangian, remains a formidable challenge despite decades of work. 
The difficulty is that the non-perturbative low-energy dynamics of strong interactions plays an 
essential role [1] , yet it has turned out to be difficult to reduce the systematic errors of lattice 
studies of these effects [2, 3] to a tolerable level, because of the prohibitive cost of treating 
chiral symmetry, final state kinematics, and physical quark masses without compromises (for 
the current status and references, see ref. [4]). 

A somewhat less challenging approach amounts to abandoning the direct computation of 
K ^ TTTT decay amplitudes in favour of determining, via lattice simulations, the low-energy 
couplings (LECs) of the effective chiral weak Hamiltonian that describes these decays [3]. 
This can be achieved by matching lattice measurements of suitable correlation functions to 
the same correlation functions computed within chiral perturbation theory (xPT). Among the 
simplifications thus achieved are that the matching does not necessitate physical kinematics, 
and that physical quark masses are not needed either, as long as the regime of validity of 
XPT is reached. This however requires sufficiently large volumes and small quark masses. 

There are various possibilities for the order in which the volume is increased and the quark 
masses are decreased. In fact, it turns out to be useful to approach the chiral limit by first 
decreasing the quark masses. The reason is that in this parameter range, referred to as 
the e-regime of [5] (see also ref. [6]), it is possible to work out next-to- leading order 

corrections in without introducing any more LECs than at the leading order, thereby 
putting the xPT side of the matching well under control. 

At the same time, carrying out lattice simulations in the e-regime is quite demanding. 
Fortunately, the advent of Ginsparg- Wilson formulations of lattice fermions [7]-[14], which 
possess an exact chiral symmetry in the limit of vanishing quark masses, and many subsequent 
developments on the numerical side [15], have made such simulations possible. 

In ref. [16] a strategy based on these methods was proposed, with the goal of revealing the 
role that the charm quark mass plays in K ^ tttt decays. The first step is the determination 
of the leading-order weak LECs in a theory with a light charm quark, that is in a four-flavour 
theory with an exact SU(4) chiral symmetry in the valence sector. The first results of this 
computation, from simulations in the e-regime, have been presented in ref. [17]. The next 
step of the strategy is to increase the charm mass and monitor the LECs as we move towards 
a theory with an SU(3) flavour symmetry [16, 18]. 

The observables used for performing the matching between lattice QCD and the chiral ef- 
fective theory in ref. [17] were three-point correlators of two left-handed currents and a weak 
operator. (The xPT side for these observables has also been worked out for SU(3) [19].) In 
this paper we pursue the same goal by means of a different type of observable. Indeed, we 
propose to consider a correlation function of two pseudoscalar densities and a weak operator. 
The peculiarity of this correlator is that it has poles in 1/m^ in the e-regime, when evalu- 
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ated in sectors of non- vanishing topological charge (which we define a la ref. [20]). As we 
will show, the residues of these poles are easier to compute numerically than the correlation 

functions themselves, since they typically require fewer quark propagators, which are substi- 
tuted by projectors on the zero-mode wave functions. We then use the residues to perform 
the matching between the fundamental and effective theories and determine the SU(4) weak 
LECs. 

Let us stress that while the weak LECs themselves are universal, and can in principle be 
determined with any method, it is difficult to know a priori which of the multitude of possible 
strategies is the optimal one in practice. This depends, on one hand, on the numerical cost 
of the measurements involved, and on the other, on how well chiral perturbation theory con- 
verges for the observable in question. We will compare the two methods mentioned (ref. [17] 
and the present one) on both accounts. The hope is that by carrying out this comparison 
in the quenched approximation, it will become clear whether or not a particular strategy is 
superior when one moves to the more expensive unquenched environment, or if it remains 
the best policy to probe the LECs by a combination of independent techniques. 

The paper is organized as follows. In Sec. 2 we introduce the observablcs to be computed 
in the fundamental theory and present the results of a next-to-leading order computation of 
the same observables in xPT, in the e-regime. Besides the three-point functions previously 
mentioned, we will consider two-point functions that we use for normalization. In Sec. 3, we 
present the results of a quenched numerical computation of these amplitudes in a 2 fm box 
and a new determination of the weak LECs in the SU(4) limit. We conclude in Sec. 4. 

2. Low-energy couplings from zero-mode wave functions 

In the e-regime and in a fixed topological sector, correlation functions involving quark propa- 
gators may contain poles in l/{mV)"-, where n is some integer, whenever the contribution of 
the zero-modes to the spectral representation of the quark propagator gives a non-vanishing 
contribution to the correlation function. For a number of reasons the residues turn out to 
be easier to compute than the correlation functions themselves. The idea, explored in detail 
in ref. [21], is then to use the residues of the topological poles to perform the matching be- 
tween QCD and xPT, instead of the full correlation function. Given a correlation function 
Cjy(xi, X2, ...), the residue can be isolated by 

Res 

C^(xi,X2,...) = , 'J' + ... , Res„ = lim (mF)"Cj.(xi,X2,...) . (2.1) 

[m V m^O 

In ref. [21] the two-point function of the pseudoscalar density was considered in this context. 
The presence of a pole in l/{mV)'^ implies that the corresponding residue can be computed 
fully in terms of the zero-mode wave functions: no propagator computation is required. On 
the effective theory side, the same pole does appear and, up to a certain order, the residue is a 
function of only the pseudoscalar decay constant, F, and the volume (in the quenched theory 
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additional couplings appear) . An exploratory numerical study in the quenched approximation 
was presented and the principal usefulness of the method to extract the low-energy coupling 
F was confirmed. Similar investigations have also been reported in refs. [22]. 

In the present work, we extend this idea to the computation of three-point functions from 
which the weak low-energy constants can be determined. In particular, we consider three- 
point functions of two pseudoscalar densities and a weak four-fermion operator. It is easy to 
see that such correlation functions do have poles in l/{mV)^ when computed in non-trivial 
topological sectors in the e-regime, as we now show. 

2.1. Correlators in the fundamental theory 

Following the strategy of refs. [16, 17], we consider a theory with four light and degenerate 
flavours such that = = = rric = m, which we refer to as the GIM limit. After inte- 
grating out the to first order in the weak coupling, g^, the resulting Weak Hamiltonian 
is given by 

9l 



where the operators Qf transform in the 84 and 20 representations of SU(4): 



Qf = {[OiUud ± [OiUdu) -(u^c), (2.3) 

[Ollrsuv = i^r7^^P-^u)i^s7^.P-^v) • (2.4) 

Here /cf are Wilson coefficients at the scale M]y, and Z^^ are the corresponding renormal- 
ization factors. We will follow the renormalization prescription of refs. [16, 17], that is we 
will use the RGI scheme, in which these factors are known non-perturbatively [23]. For any 
unexplained details we refer the reader to these references. 

Deep in the non-perturbative regime this effective Hamiltonian admits an expansion in 
terms of the Goldstone boson fields and can be represented as 

= ^KsVud {gt Qt + 9i Qi + •••} , (2-5) 

where are operators made out of the Goldstone field, and terms of higher order in the 
chiral expansion have been omitted. Our task is to match for the coefficients gf in the chiral 
limit, by comparing lattice simulations with xPT predictions. 

Now, given that QCD dynamics itself respects chiral symmetry, the results of such a match- 
ing are independent of the precise flavour indices appearing in eq. (2.3), as long as the op- 
erators remain traceless and have the correct symmetry properties. In practice, it is indeed 
convenient to consider the operators 

of = [Oijrsuv ± [OljrsTO , (2.6) 
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with all indices different, rather than Qf ; these operators are automatically traceless, and 
no subtraction of the type in eq. (2.3) is needed. The results for the matching of gf are 
nevertheless guaranteed to be identical. 

We want to carry out the matching in a finite volume through the computation of the 
following bare three-point functions in a fixed topological sector of charge v. 

Atixo - zo, yo -zo) = - lun imV)^ [ [ {d,,P\x)Of {z)dy,P\y)) , , (2.7) 

where the bare pseudoscalar density reads P° = i^^^T""^, and m is the bare quark mass. 
These amplitudes get contributions from two possible contractions, colour-disconnected, Ai^, 
and colour-connected, A^: 

Ai: (TO. rph I rpa rph I rpa rpb I /TiO nib \ A -1- A (O &\ 

■^v — K-'-ur-'-vs'^ -'-vs-'-ur^ -'-vr-'-us^ -'-us-'-vr) -f^v ^ -f^v , {'^■o) 

where 

A^ = limimVf / d^ody^Ti [Sm{x, z)'yi^P-Smiz, x)'y5]Ti [Sm{y, z)'y;^P-Sm{z,y)'y5] , 

^x Jy 

A^ = -\\m.{mV)^ / / 5a;o9j^oTr [5m(a;, 2:)7^P_S'm(2;,?/)75S'm(?/,z)7/,P_5;„(2;,x)75] . 

7x Jy 

(2.9) 

Here Sm is the massive quark propagator. 

It is convenient to normalize these three-point functions with bare two-point functions of 
the form 

-i'IV[T"T'']5,(xo-^o) = Xyai{mV) I {d^^P\x)L\{z))^ , (2.10) 
where the bare left current reads Lq = 5'7oP_r"^. Carrying out the contractions, we get 

B^{xo-zo) = lim{mV) [ d^,Tr[Sm{x,z)^oP-Sm{z,x)j5]- (2.11) 

Note that the two-point function in eq. (2.10) can be related through the non-singlet axial 
Ward identity to the two-point function of two pseudoscalar densities, considered in ref. [21]: 

TV [T"T^]ZaS,(xo -zo) = i lim {mV)d,, [ {P''{x)ZaLI{z)), 

= -i \im {mV)d,, I {P%x)ZaLI{z)), 

"^-^0 Jx 

= lim (mV) / (P"(x) P\z))^ , (2.12) 

where Za denotes the renormalization constant of the currents Lq. Here we made use of the 
fact that the product mP°' does not require renormalization. 
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Now, let us see why it is useful to consider the quantities in eqs. (2.9), (2.11). The point 
is that Ai, , Ai, and Bi, are zero, unless some of the quark propagators are saturated by zero 
modes. Let us denote by Vi{x) € fC the zero- mode wave functions (recall that the dimension 
of the kernel of the Dirac operator is dim(/C) = that are normalized as 

/ vl{x)v^{x) = V . (2.13) 

Jx 

The spectral representation of the quark propagator then reads 

s„(x,rt=x;^5M^lM + .... (2.14) 

i=l 

We can define the sources 

r]i{z;xo) = P-y-Smiz,x)P^Vi{x) , 

Jx 

vjiz-xo) = -d^, [ vl{x)P^Sm{x,z)P-^, (2.15) 

where x is the chirality of the zero-modes. Given that ^y^P- = P+^y^P-, we note that, 
depending on chirality, only some of the propagators can be saturated with zero-modes: for 
u < 0, the ones multiplying P_, and for u > 0, the ones multiplying P+. If > 0, the 
disconnected and connected amplitudes of the three-point functions thus become 



Ai,{xo - zo, yo - zo) = Hin / vl{z)^^r)i{z; xq) ^ v]{z)jij,rij{z; yo))^ , 

""^ "'^ i=i j=i 

1 f I'^l 

A^{xo- zo,yo- zo) = -Jim J (^"^ vliz)'y^rij{z;yo)vj{z)-ff,rii{z;xo)'^ ^ , (2.16) 



while for < we get 

A,,{xo- zo,yo- Zo) = hin / Q2r]l(z;xo)jfj,Vi{z)^r]j{z;yo)jfj,Vj{z 

•''^ i=i j=i 

A^{xo- zo,yo- Zo) = -hin -3 U^ri\{z\xo)lij,Vj{z)ri]{z\yo)lij,Vi{z)'^ -(2.17) 

For the two-point function of eq. (2.11), the positive chirality case v > Q yields 

1 f '''' 

B„{xo - Zo) = Jim J (^^vl{z)'yor]i{z;xo)'^^ , (2.18) 



5 



while for < we arrive at 



Bt,{xo - zo) = lim — 3 




^rjl{z;xo)joVi{z) 



i=l 



V 



(2.19) 



Finally, the Ward identity of eq. (2.12) implies 




(2.20) 



Here the limit m — >^ has been taken analytically on the right-hand side, while it needs to be 
taken numerically on the left-hand side (cf. eqs. (2.18), (2.19)); therefore eq. (2.20) provides 
a non-trivial test on our ability to approach the limit needed in cqs. (2.16)-(2.19). 

It is clear from eqs. (2.16)-(2.19) that a number of inversions equal to twice the topological 
charge, i.e. 1\v\ (since xq and yo need to be fixed), is sufficient for constructing the correlation 
functions, whilst averaging over all the spatial positions of the three sources. When employing 
the method of refs. [16, 17], which is based on the left-handed current, the summation over 
the spatial positions of the three sources is only possible through low-mode averaging (LMA) 
[24, 25], and only for the contributions of the low modes. The price of LMA is 12 -|- 2 x A^iow 
inversions, where Mow was the number of low modes treated separately. Typically Mow can 
be as large as 20 for 2 fm boxes, and hence the numerical cost can be quite substantial. 

On the other hand, if the low modes of the Dirac operator are known, as for example would 
be the case if low-mode preconditioning [15] is used, it is possible to perform an additional 
averaging over time translations for the low-mode contributions to the correlation functions 
defined above. It is important to stress however that this extra low-mode averaging does not 
involve any additional inversion, therefore the overhead is not proportional to Nioyj, as in the 
standard case [24]. We will describe briefly how this works in section 3.1. 

Summarizing, in order to perform the matching between the fundamental weak Hamilto- 
nian and the effective one, we will be considering the bare ratios 



The renormalized ratios needed for matching the LECs in eq. (2.5) are then obtained by 
multiplying these correlators by the renormalization factors Z^i/Z\ (cf. eq. (2.41) below); 
the procedure is identical to the one in ref. [17] and we refer to that reference for details. 

2.2. Correlators in chiral perturbation theory 

We now present the results for the observables just deflned in the Chiral Effective Theory. 




At,{xo - zq, yo - zq) ± A^{xo - zo,yo - zp) 

B^{xo - ZQ)B^{yQ - Zq) 



(2.21) 
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2.2.1. Two-point function 



The two-point correlation function that needs to be computed in xPT, corresponding to 
eq. (2.10), is given by 

-iTr [T"T^]H,(xo - zo) = lhn (mV)d,, [ (P«(x) J^iz)), , (2.22) 



where 



rra 



— Tr 
2 

i-Tr 
2 



T'^Udi.U^ 



+ ... , 
+ .. 



(2.23) 
(2.24) 



and contact terms have been omitted from eq. (2.22). 

Since the mass is taken to zero in eq. (2.22), the computation is carried out according to 
the rules of the e-expansion [5]. We work up to next-to-leading order (NLO). The results for 
the contributions of the individual graphs, as well as the various zero-mode and spacetime 
integrals appearing, are listed in appendix A. 

Defining = (xq — zq)/T and p = T/L, and considering the unquenched theory, the final 
result from eq. (A.29) becomes, after replacing E{x) = G{x)/Nf and using eq. (A.27), 

TB,{xo-zo) = IH{l + ^^(lH + ^)^i(r.)} , (2.25) 



where 



hiir) ^ 



(r mod 1) 



I 



H 2 



(2.26) 



In Fig. 1 we show this result for different values of Iz^l in a symmetric box of size T = L = 2 fm. 
The LO results are constants at so that all time dependence results from the subleading 
corrections. 

In the quenched case we consider the e-counting described in rcf. [21], introduced in order 
to set up a paramctrically convergent perturbative series. As explained there, three new 
couphngs can in principle appear at NLO: a,mQ,if. The effect of a is to replace 1/Ni — > 
a/2Nc; however, this contribution is suppressed in the counting of ref. [21] and will be omitted. 
Similarly, contributions from mg are also suppressed and omitted. On the other hand, effects 
from the coupling K could be of order unity, and need to be studied explicitly. 

The coefficient K has two main effects. First of all, the pseudoscalar density is modified to 



pa _ 



i—Ti 
2 



T^iU - U^) 



K$oTr 



T'^iU + U^) 



+ ... , 



(2.27) 



where we have for brevity kept the notation of the unquenched theory, usable in the replica 
formulation of quenched [26]; for the notation in the supersymmetric formulation [27], 
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Figure 1: The two-point correlation function TBy{xQ — zq) for = 2 (solid) and A^f = 
(dotted) as a function of = {xq — zo)/T, for T = L = 2 fm. The pion decay constant F 
has been fixed to 93 MeV for Ni = 2 and to 110 MeV for iVf = 0. 



see eq. (3.2) of ref. [21]. In eq. (2.27), ^>o = fTr[-ilnC/] is the singlet field. Second, the 
zero-mode integration measure is modified: 



/uoeU(N)(-) def'f^o exp 


"fTr([/o + C/J) + 2i.l^TV(C/o 


- U^o)_ 


/uoeU(N) def^f/o exp 


2 


TT{Uo + U^o) + 2u^Tt{Uo- 







= ^ 77^ TT^ ■ (2-28) 



To first order in K, the effects come from the standard tree-level contribution computed with 
the modified weight of eq. (2.28), and from a new tree-level term containing the X-term from 
eq. (2.27). 

The results for the two new contributions are given in eqs. (A. 14), (A. 15) of appendix A. 
However, inserting the zero-mode integrals from eqs. (A. 24), (A. 25), these contributions 
cancel against each other. Therefore, 

TBl{xo-zo) = |j.||l + ^^/,i(r,)| , (2.29) 

and the result is very close to that in the full theory (were it not that F is different). Some 
examples are shown in Fig. 1. 
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2.2.2. Three-point function 

Next we consider the three-point function corresponding to eq. (2.7) in : 



At{xo - ZQ,yo- zq) = - lim (mV)' 

m—*0 



Here, hke in eq. (2.6), 

where the weak operators are given by 

^4 



Or 



= — (d^UU^) (d^UU^) + ... . 

4 \ / ur \ / vs 



(2.30) 
(2.31) 

(2.32) 



Again, contact terms have been omitted from eq. (2.30). 

Like for the two-point function, we work up to NLO. The results for the contributions of 
the individual graphs, as well as the various zero- mode and spacetime integrals appearing, 
are listed in appendix B. 

The result obtained after summing together all the graphs can be written as 

Aziz 'T^ 4_ 7^^ -1- 7^*^ -I- 7^*^ 7^^ 

"^u \ ur vs ' vs ur vr us us vr 



where, inserting E{x) = G{x)/Nf into eq. (B.33) and using eq. (B.28), one obtains 
Ay±Ay= (^1 =F j^) |.Sj.(xo - zo)By{yo - zq) 
lb 



(2.33) 



2/3i/9-5 + /i(r^) + hiTy) - /ii(r^) - /ii(r^) + ( 1 T —)H{T^,Ty) \ . (2.34) 



Here Tx = {xq — zo)/T, Ty = (yo ~ zo)/T, p = T/L, (3i is a shape coefficient depending on the 
value of p [28, 29], the function hi is defined in eq. (2.26), and 

H{Tx,Ty) = h[{Tx)h[{Ty) - hi{Tx - Ty) - h[{Tx) - h[{Ty) h'i{Tx - Ty) , (2.35) 

fi{r) ^ [h[{T)]' + Y,[\pfCp{r)' + C^{r)'] . (2.36) 

Furthermore, p = 27rpn with n = {ni,n2,'n^) being a vector of integers, and 

^ / ^ _ cosh{|p| [(rniodl) - i] } _ sinh{|p| [(rmodl)-^]} 

^''^^>= 2|p|sinh(|p|/2) ' ^^^^>- 2sinh(|p|/2) " ^^"^^^ 

The first term in eq. (2.34) has the form of a factorized contribution. We can cancel out 
this term by taking the ratio of the three-point function with respect to the product of two 
two-point functions, in analogy with eq. (2.21): 

Au{xo - zo, yo - Zq) ± Au{xo - zo,yo - zq) 



By{xo - zo)By{yo - zq) 



= (iT p)[l±r±(zo) 



(2.38) 
(2.39) 
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We have adopted a notation here where xq, yo are assumed fixed, so that r± is a function of 
zq only. Inserting the tree- level result Bi, = \u\/T in the numerators when dividing the NLO 
correction, we then obtain from eq. (2.34) that 

r±(zo) = [2/3ip-i + /i(t,) + /i(t,) - /ii(r,) - /ii(r,) + (l ^ ]^)^(^- ^v)] ■ (2-40) 



Note that r± is independent of u, i.e. the topology and volume dependences have completely 
factorized at this order. The low-energy couplings gf can now be obtained from the matching 



gtnt = kt^Rt, (2.41) 

at sufficiently large distances between the sources. 

The quenched result in the counting of ref. [21] is obtained by simply leaving out the term 
l/7Vf, because terms involving a and ttZq are of higher order, like for the two-point correlator, 
and effects from the coupling K cancel at this order, as demonstrated in appendix B (cf. 
eqs. (B.13), (B.14), (B.20), (B.21)): 



rlizo) 



{FLf L 



Wlp'^ + h{Tx) + h{ry)-hi{T:,)-hi{Ty)+H{T:,,Ty) . (2.42) 



The most efficient way of evaluating numerically the amplitudes and is by fixing the 
temporal position of the sources, xq and uq, so that the three-point correlator is measured 
as a function of the temporal position of the weak operator, zq. In order to maximize the 
separation between the three operators we take xq ~ T/3 and uq ~ T — xq. The signal should 
be best when the weak operator is near the origin, and simultaneously the NLO correction 
r±{zQ) is minimized. The corresponding values of (FL)^rj.(0) for various box shapes are given 
in Table 1. There is a very strong dependence on zq (cf. Fig. 2) which results from the fact that 
the functions fi{Tx) and fi{Ty) arc divergent at r^- = and Ty = 0, respectively, while they 
fall off exponentially, ~ ex];){—2TTpTx), ~ exp(— 27r/9Ty), away from these points^. Thereby 
a way to decrease these corrections is to increase p = T / L . In this respect the situation 
is opposite to that in ref. [17] where left-handed currents appear in place of pseudoscalar 
densities; then the NLO corrections increase rapidly with r/L>2 [16]. 

Fig. 2 shows the result for 1 it for various box volumes, for xq = T/3, yo = 2T/3. 
Unfortunately. NLO corrections seem to be rather large at L = 2 fm. 



^These terms could be decreased by increasing r^, Ty towards 1/2, since nothing dramatic happens as 
Tx — > Ty according to NLO expressions (cf. Table 1). However, this situation may be specific to the SU(4) 
case, where there are no quark propagators connecting the pseudoscalar densities, while in the physical SU(3) 
case it is probably reasonable to keep the pseudoscalar densities somewhat apart from other. 

^Note however that increasing p takes us closer to the so-called J-regime [30]. 
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V 


xo/a 


yo/a 






(FL)Vi(O) 


16* 


4 


12 


0.25000 


0.75000 


0.89309 


16* 


5 


11 


0.31250 


0.68750 


0.62839 


16* 


6 


10 


0.37500 


0.62500 


0.51281 


16* 


7 


9 


0.43750 


0.56250 


0.41875 


16^32 


11 


21 


0.34375 


0.65625 


0.45251 


16^32 


12 


20 


0.37500 


0.62500 


0.43174 


16332 


13 


19 


0.40625 


0.59375 


0.38822 


16^32 


14 


18 


0.43750 


0.56250 


0.32157 


16^32 


15 


17 


0.46875 


0.53125 


0.23162 


24* 


8 


16 


0.33333 


0.66667 


0.58270 


24* 


9 


15 


0.37500 


0.62500 


0.51281 


24* 


10 


14 


0.41667 


0.58333 


0.45099 


24* 


11 


13 


0.45833 


0.54167 


0.38435 



Tabic 1: Examples of values of (FL)^r!^(0) for various box shapes. Note that in the quenched 
hmit, r^(0) is independent of the channel ±, cf. eq. (2.42). 





/3 


ro/a 


V 


A^low 


|z/| 


conr 


a ' a 


Za 


am 


Ai 


5.8458 


4.026 


16* 


20 


1-5 


118,94 

99,73,65 


5,11 


1.710 


0.0015,0.0025, 
0.005 


A2 


6.0735 


6.072 


24* 


20 


2-5 


92,63 
51,55 


8,16 


1.505 


0.002, 0.0033, 
0.0067 



Table 2: Parameters of the simulations. 
3. Numerical results 

We have carried out a numerical test of the method outlined above, in the quenched approx- 
imation. Table 2 shows the simulation parameters. We have considered a symmetric lattice 
T = L ~ 2 fm at two different lattice spacings in order to test for scaling violations in these 
observables. 

To evaluate the correlation functions of eqs. (2.16)-(2.19), we have computed chiral quark 
propagators on quenched background gauge configurations, using the Neubcrgcr-Dirac op- 
erator with s = 0.4. For all details of the numerical implementation we refer the reader to 
refs. [15, 24, 16], whose techniques we adopt. 

Before presenting the actual data, let us describe how low-mode averaging (LMA) can be 
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Nj, = 2 Nj = 




Figure 2: The functions lzizr± (eq. (2.39)) for Nf = 2 (left) and A^f = (right) as a function 
of To = zq/T, for xo = T/3, yo = 2T/3. The pion decay constant F has been fixed to 93 MeV 
for iVf = 2 and to 110 MeV for A^f = 0. The sohd hne corresponds to T = L =2 fm, the 
dashed to T = L = 3 fm, and the dotted to T/2 = L =2 fm. 



used to reduce statistical fluctuations in the signal. 
3.1. Low-mode averaging 

Our observables, eqs. (2.16)-(2.19), involve both zero-mode wave functions and actual quark 
propagators. Since we employ low-mode preconditioning [15] for determining the quark prop- 
agators, which requires the computation of a few lowest eigenvectors of the Dirac operator, we 
can use these in order to perform an extra averaging over time translations of the low-mode 
contribution to eqs. (2.16)-(2.19). In other words, the LMA technique [24, 25] gives only an 
extra averaging over time in our case, but is nevertheless helpful as we will see, since the 
numerical overhead involved is negligible. 

The main idea of LMA is to substitute the chiral propagator in eqs. (2.16)-(2.19) by 

P.^SUz,x)P^ = ^ vI/,(z)®^^(x) + P_^5^"^z,x)P^ , (3.1) 
fc=i 

where S^^{z,x) is the inverse of the massive Dirac operator in the subspace orthogonal to 
the eigenspace of the approximate low-modes, and the chiral components of ai"e given by 

P-^^k{x) = Wk{x), P^^kix) = —P^j5DP-^Wk{x) , (3.2) 
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where Wk(x) are the approximate eigenfunctions of the operator D^D for the eigenvalue A^, 
while Ofc = \/Ajfc. We indicate how this works with the two-point correlator. 
The LMA evaluation of eq. (2.18) is based on the separation 



where the "high-mode part" reads 



Bu — bI + B^, 



1^1 



with 



B^{xo-zq) = hin j (^vl{z)-iQf}i{z-xo)^ , 

-^^ 1=1 

while the low-mode contribution is 

W\ ^low . . 

Bi{t) = lim VEv/ 5{xo-zo-t)lvl{z)^oP--^k{z)d,,bl{x)P+v,{x) 
i=i k=i -^^'^ 



(3.3) 



(3.4) 



(3.5) 



(3.6) 



The LMA evaluation of the three-point function is carried out in complete analogy. 



3.2. Two-point function 

In Fig. 3 we show results for the two-point correlator Bi,{t) (cf. eqs. (2.18), (2.19)) in different 

topological sectors at the lightest quark mass. The open/full symbols corresponds to the 
results without/with LMA. There is a strong dependence on |z^|, as expected from xPT. 

Since there is clear evidence for NLO corrections, we consider a two-parameter fit of the 
form 

TB,{t) = a, + 2pM (r) , r = | , (3.7) 

where hi is from eq. (2.26). The temporal dependence in all sectors is perfectly compatible 
with the function /ii(r) as illustrated by the solid lines in Fig. 3, which are the results of the 
fits in the time interval = 5a — 11a. The small dependence on the quark mass is perfectly 
linear, so the vahics of a^, and f3i, are linearly extrapolated to the zero mass limit. The results 
for ay and [3y in the chiral limit are summarized in Table 3. The jackknife procedure has 
been used to estimate the errors. 

In Fig. 4 we show the results for and (3y as a function of together with the NLO 
expectations. In the case of a^,, the NLO prediction = \v\ is extremely well reproduced at 
the per cent level. The prediction for (3i,-, on the other hand, depends on F. The dashed line in 
the figure corresponds to a fit to the NLO prediction leaving F as a free parameter. The best 
fit values are (FL)ai = 1.19(2) and {FL)a^ = 1.14(2) with x^/d.o.f- 9 and 12, respectively. 
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..8 ^^^~tr-tr-&^ 8--v=2 
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t/a 



Figure 3: The measured TBy{t) for lattice Ai and quark mass am = 0.0015. Open/full 
symbols correspond to data without/with LMA. Error bars are in most cases smaller than 
the symbol sizes. 



which are rather bad. Clearly the dependence is not properly reproduced at NLO, however 
it seems that the discrepancy could be ascribed to higher order chiral corrections. 

We have seen that the Ward identity, cq. (2.12), relates Bu(t) to the topological zero-mode 
contribution in the correlator of two pseudoscalar densities. This quantity was studied in 
ref. [21], and actually the chiral corrections were computed to one order higher than here. 
Although the expressions are rather complicated and involve new time-dependent functions, 
a convenient way to try to include these corrections is to consider a Taylor expansion around 
the middle point, t = T/2. Prom the results of ref. [21], we expect that 

ri?,(i)=7, + <^. +... , (3.8) 

where^ 

_A,-|_i_(^+..^^2(.^))+|]}. (3.9) 

For a symmetric box with T = L the shape coefficients read p = 1, /?i = 0.14046098, and 
71 = -0.05712765. 

Numerically the results for 5i, in a fit of the form in eq. (3.8) arc identical to those for /3jy 
in Table 3. The solid line in the lower plot of Fig. 4 corresponds to the prediction of eq. (3.9) 

^We have set a = 0, since this parameter was consistent with zero in the analysis of ref. [21]. 
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(FLf 
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\u\ a,(Ai) (Ai) a,(A2) /?. (A2) 



1 


1.00(1) 


1.4(3) 






2 


1.98(2) 


3.9(4) 


2.02(2) 


4.5(4) 


3 


3.02(3) 


7.3(4) 


3.00(3) 


8.8(5) 


4 


4.00(4) 


11.6(5) 


3.97(4) 


11.9(5) 


5 


5.01(5) 


15.8(6) 


5.00(4) 


17.7(6) 



Table 3: Results for and (3^ from a fit to eq. (3.7) of the measured two-point correlator, 
subsequently extrapolated to the zero-mass limit. 

using {y"^) = = 0.059y/rQ [31] and leaving i*" as a free parameter. The result of the fit 
is (FL)ai = 1.12(2) and {FL)^^ = 1.07(2) with xVd-o.f = 0.2, 1.3, respectively. The resuhs 
for F are essentially the same as the central values in ref. [21] (this is almost trivial in the 
light of Sec. 3.3 below)^; however, the current error bars are much smaller, because higher 
topological sectors which are less noisy have been considered. Besides, we have not assigned 
any errors to a and Xt, as wc did in [21]. The results are in reasonable agreement also with 
the determination of F from the left-current two-point correlator in ref. [24] , although it must 
be noted that the box size was significantly smaller there. 

3.3. Ward identity and chiral extrapolation 

The Ward identity of eq. (2.12) is a good test of the extrapolation m ^ needed for the 
correlator B^(t). 

In Fig. 5 we show results for the ratio Za = Di,{t)/Bi,{t) (notation as in eq. (2.20)) as 
a function of am for the different topological sectors, normalized to the value Za obtained 
by conventional means in ref. [32]. In the limit m ^ 0, the ratio should approach unity in 
all topological sectors. The level of agreement between the different sectors and with Za is 
shown in Table 4. Given that only the zero-mode contributions to both sides of the Ward 
Identity are included, this is a strong check of the whole procedure, and furthermore indicates 
that the small residual extrapolation to zero quark mass is under good control. 

3.4. Three-point function 

In Fig. 6 wc show one example of a Monte Carlo history for the three-point functions Ay^A^,, 
eqs. (2.16), (2.17), with and without LMA. Clearly LMA improves the signal significantly. 

''Note that a similar physical box size was used in ref. [21] as in the present study. 
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V 



Figure 4: Top: ayl\v\ versus \v\ for lattice Ai (left) and lattice A2 (right). The dashed line 
is the NLO expectation. Bottom: fiy versus \v\ for the same lattices. The dashed line is 
the best fit NLO prediction, the solid line is the best fit NNLO prediction. In both cases 
open/full symbols correspond to without/with LMA. 







Za/^a(A2) 


1 


1.000(6) 




2 


1.009(4) 


0.991(3) 


3 


0.999(4) 


0.995(3) 


4 


1.008(3) 


0.999(3) 


5 


1.000(3) 


1.002(3) 



Table 4: The ratio Za = D^jB^ (cf. eq. (2.20)), obtained from the saturation with zero modes 

of the Ward Identity (with LMA treatment of the non-zero mode part of By\ normalized to 
the conventionally determined Za [32]. The errors do not include the error on Za which is 
about 3 per mille. 

The improvement is more pronounced for the smaller topologies and masses as expected. In 
all of the following we consider only the LMA results. 

In Figs. 7, 8we show the results for the bare ratios i?^(xo — 2;0)yo — ^o)/(lTl/|'^|) on lattice 
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Figure 5: Chiral extrapolation of the time-averaged ratio D,^(t) / Bi,{t) (cf. eq. (2.20)) nor- 
malized to Za from ref. [32] , for lattice Ai . 



Ai as a function of r = zq/T, at fixed Xq = 5a, yo = 11a. The quark mass is am = 0.0015. 
There is a clear signal near r = 0. However, the temporal dependence docs not seem to be 
as pronounced as expected from NLO xPT. The pattern is similar for the lattice A2. 

An interesting combination to consider is the product R^R~ , since NLO corrections cancel 
in this quantity (cf. eq. (2.39)). Writing the weak LECs as 

9t = [gtr-^ ^ , (3-10) 

eqs. (2.39) and (2.41) imply that we may expect: 

R-^Ru = [gtgi]'''"' (i - ^) + • • • ■ (3-ii) 

The results for the lattice Ai are shown in Fig. 9. A constant fit around zq = 0, followed 
by a chiral extrapolation in each topological sector, gives the values shown in Table 5. The 
agreement with the result obtained from the left-current three-point functions [17] is quite 
good. Note that the renormalization factors for the lattice Ai arc the same as those in 
ref. [17] and therefore the bare couplings can be compared directly. Those for the lattice A2 
have not been evaluated. However, the difference is expected to be well below the statistical 
uncertainty. Indeed a LO computation in bare perturbation theory yields variations at the 
1-2% level, and the results of ref. [23] indicate that this perturbative estimate is realistic. 

In order to then estimate [^j'']^^''^ and [^f J^^"^^ individually, we consider the following fits 
of the LMA data: 
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25 50 75 100 



Figure 6: The Monte Carlo history of the three-point correlators vlj/ — Ai, (top) and A^^ + A^ 
(bottom), at zq = 0, normalized to the average, on the lattice Ai, for = 3 and am = 0.0015. 
The dark hne corresponds to the LMA and the light one to the non-LMA amplitude. 



IH 






[gigi]""^'' [17] 


2 


0.73(53) 


0.72(39) 




3 


0.94(20) 


1.10(34) 




4 


1.37(20) 


1.60(36) 




5 


1.64(20) 


1.50(17) 




w.a. {\v\ > 2) 


1.32(11) 


1.45(14) 


1.47(12) 


xVd.oi. 


3.1 


0.7 





Table 5: Values of the product of the bare couplings obtained from fits to a constant in the 
time interval j^oj < a for lattice Ai and jzoj < 2a for A2. The last row shows the weighted 
averages (w.a) over topological sectors jz^j > 2 for the two lattices and also the result obtained 
in ref. [17] with a different method. 
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Figure 7: R^/{1 — in different topological sectors for the lattice Ai as a function of 

zq/T for fixed xq = 5a, yo = 11a, at the smallest quark mass am = 0.0015. The horizontal 
lines represent the la boundaries of a LO fit, while the curved line is the best NLO fit. 



Fit A. At LO we expect (cf. eq. (2.39)) 



It 



(3.12) 



therefore we fit / ^1 =p j^j to a constant around zq = 0. The results are shown as the 
bands in Figs. 7, 8. 

The results for the bare couplings are listed in Table 6. The results on the lattices Ai and 
A2 are again perfectly compatible, which implies that scaling violations are well below the 
statistical uncertaintities. There is a significant difference with the [fi']^]'^'''^'' obtained from a 
LO matching of left-current three-point functions in ref. [17]: while our [fi'j"]'^'*'^'' is smaller, 
our [s'J^]^^'^*^ is larger. This is, however, not unexpected, given that chiral corrections tend to 
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Figure 8: R~ /{I + l/|z^|) in different topological sectors for the lattice Ai as a function of 
zq/T for fixed xq = 5a, yo = 11a, at the smallest quark mass am = 0.0015. The horizontal 
lines represent the la boundaries of a LO fit, while the curved line is the best NLO fit. 



decrease [g^ J'^^''*' and increase [g'j'']'^'^'^'^ in ref. [17], and do the opposite in the present work. 
Fit B. At NLO we expect (cf. eq. (2.39)) 



[at 



ztnbare 



1^ M 
\l^\ 



l±4(zo) 



(3.13) 



Numerical values of (FL)^r^(2o) for small zq are shown in Table 7. Taking a value FL ~ 1.1, 
the NLO prediction gives a stronger temporal dependence than that seen in the data. The 
fits have a bad if more that 3(5) points are included for lattice Ai(A2). This already 
indicates that NNLO could be significant. With so few points it does not make sense to do 
a two parameter fit, leaving the normalization of the NLO correction free, because there is 
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Figure 9: R^R~/{1 — l/i^^) in different topological sectors on the lattice Ai as a function 
of zq/T for fixed xq = 5a, yo = 11a, at the smallest quark mass am = 0.0015. The bands 
represent the la fits to a constant. 



not much curvature in the data. Therefore we fix (FL)'^ = 1.14 — 1.40, corresponding to the 
range of values obtained from fitting various two-point functions here and in ref. [33]. We 
then perform linear one-parameter fits for [^j'^]'''^'^'^. The results are summarized in Table 8. 

The results obtained for [s'j'^]'^'^'^'' on the lattices Ai and A2 are again fully compatible, 
indicating small scaling violations. However, the effect of a 10% uncertainty in F results in 
a significant systematic uncertainty especially in [gf]*''*'^'', which is much more sensitive to 
NLO corrections. In fact the difference between the results for [^^f ]'^'^'^'' obtained from LO 
and NLO matchings is too large for even the latter results to be trustworthy. 

A comparison with the results of a NLO matching in ref. [17] shows that [5]^]'^'*'^° is in rather 
good agreement. This is quite non-trivial, given the very different NLO chiral corrections in 
the two cases, and could be an indication that NNLO corrections are not very important for 
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[g-i 



ibcire 



[9t 



+lbare 



1 1 


Ai 


A2 


[171 

L J 


Ai 


A2 


[171 

L J 


1 


3.1(1.6) 












2 


0.8(7) 


1.7(0.7) 




0.9(2) 


0.44(13) 




3 


1.2(2) 


1.46(35) 




0.76(8) 


0.80(10) 




4 


1.65(15) 


1.85(31) 




0.83(6) 


0.89(10) 




5 


1.73(12) 


1.74(15) 




0.95(5) 


0.86(6) 




w.a. {\u\ > 2) 


1.60(8) 


1.72(13) 


2.42(13) 


0.87(3) 


0.85(5) 


0.60(4) 


xVd.o.f. 


2.6 


0.4 




2.4 


0.2 





Table 6: Values of the bare couplings obtained from a LO fit in the time interval \zo\ < a on 
lattice Ai and \zo\ < 2a on lattice A2. 



zo 




[(FL)Vi(^o)]A, 





0.62839 


0.58270 


a 


0.73420 


0.61454 


2a 


1.23214 


0.72759 



Table 7: Values of {FL)'^r^{zo) for zq = 0, a, 2a, with xq 
xq = 8a, Do = 16a for lattice A2. 



5a, yo = 11a for lattice Ai and 



this quantity. On the other hand, our [^f J^^"^*^ lies significantly above the result of rcf. [17], 
but this could be accounted for by an effect of 30% in the NNLO corrections, which docs not 
appear unreasonable, given that the size of the NLO corrections is 50 - 60%. It is interesting 
to note, however, that due to the fact that the there is a cancellation between the LO and 
NLO terms in R~ , the uncertainty in the NLO corrections has a bigger relative impact in 
the determination of [gf]^^"^^ than in [5]'']'^'^'''^. 

Obviously it is necessary to bring this systematic error under control, which can only be 
achieved by going to larger volumes and performing a detailed finite-size scaling study. Given 
that this is a quenched exploratory study and that going to larger volumes in the quenched 
approximation is no guarantee of success, we will not pursue this further here. We will 
consider, however, alternative estimates of [^jf J*^^^ that could be less affected by higher order 



22 



[17] 



2 


2(2)- 


1.6(1.4) 


4(2) - 3(1) 




3 


3.1(5) 


2.5(4) 


3.4(8) 2.7(7) 




4 


4.2(4) 


- 3.3(3) 


4.3(7) - 3.5(6) 




5 


4.4(3) 


- 3.4(2) 


4.0(4) - 3.2(3) 




w.a. > 2) 


4.1(2) 


-3.2(2) 


4.0(3) - 3.2(3) 2.33(11) 


xVd.o.f. 


2.5 


- 2.0 


0.4 - 0.4 












[^+]ba.e 






A 


1 


A2 


[17] 


2 


0.56(13) - 


0.60(14) 


0.30(9) - 0.32(10) 




3 


0.47(5) - 


0.51(5) 


0.53(7) - 0.57(7) 




4 


0.52(4) - 


0.55(4) 


0.59(7) - 0.63(7) 




5 


0.59(3) - 


0.63(3) 


0.55(4) - 0.59(4) 




w.a. {\v\ > 2) 


0.55(2) - 


■ 0.58(2) 


0.57(4) - 0.61(4) 


0.63(4) 


xVd.o.f. 


2.4- 


2.6 


0.3 - 0.2 





Table 8: Values of the bare couplings obtained from a NLO fit in the interval \zq\ < a. The 
ranges indicated correspond to considering {FL)"^ = 1.14 — 1.40. 

corrections. 

One possible strategy^ is to obtain [^f ]^^''^ indirectly from [g^g^^^^^ and [5]'"]'^^^. The 
first quantity is extracted from the product R^R~ where the NLO correction vanishes, while 
the second quantity is extracted from where the NLO contribution has the same sign as 
the LO, and therefore the uncertainty due to higher order effects is less relevant on relative 
terms. The results of such an approach are summarized in Table 9. The results for [yf ]^'^'^^ 
are now very close to those in ref. [17], and are significantly less sensitive to the uncertainty 
in F, as expected. 

Alternatively, as is clear from Table 1, on could decrease the chiral corrections very signifi- 
®This was the strategy followed in ref. [17]. 
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[17] 



4 



5 



2 



3 



1.3(1.1) - 1.2(1.0) 2.4(1.1) - 2.2(1.0) 

1.99(32) 1.85(30) 2.1(5) - 1.9(5) 

2.66(24) - 2.47(22) 2.71(45)- 2.54(42) 

2.79(19) - 2.60(18) 2.65(24) - 2.48(22) 



w.a. > 2) 2.61(14) 2.42(13) 2.58(20) - 2.42(18) 2.33(11) 
xVd.o.f. 2.3 - 2.3 0.5 - 0.6 



Tabic 9: Value of the bare coupling [g^ ]^'^'^^ obtained from a NLO fit of R^Rj^ and in 
the interval jzoj < a for lattices Ai and A2. The ranges indicated correspond to considering 
{FLf = 1.14-1.40. 

cantly by increasing the distances Tx and Ty between the pseudoscalar densities and the weak 
operator, at the expense of decreasing \tx — Ty\. Note that chiral corrections (up to NLO) are 
in fact insensitive to the last separation. In principle, the distance between the pseudoscalar 
sources has to be large enough compared with the cutoff and the physical distance scales of 
QCD. However it is an empirical observation that two-point functions approach the asymp- 
totic behaviour very fast in the e-regime, and we are therefore confident that it makes sense 
to investigate the three-point functions for larger values of and Ty, even on a lattice as 
small as r = 16a. Provided the effects of higher scales can be neglected, choosing xq = 7a 
and yo = 9a on the lattice Ai can reduce the chiral corrections by 30%. The results obtained 
for this choice are summarized in Table 10. 

A few observations are in order. The changes in [5^]'^^'^'^ and [gfgi]^'^^^ with respect to 
the case Tx — Ty c::i \tx — Ty\ are quite small, but the change in [yf ]'^'^'^'^ is very significant (yet 
still at the level expected from NNLO chiral corrections), bringing the value of [yf j^^"^® to 
agreement with that from the indirect determination, and with that of ref. [17]. The x^/d.o.f. 
of the fits get improved and the effect of the uncertainty in F on [5^ ]'''^'^^ gets reduced to 
the level of 10%. In general the agreement of these results with those of ref. [17] is quite 
remarkable. Unfortunately the distance between the sources is too small to be confident that 
the effect of higher scales is negligible, but these results provide further evidence that the 
discrepancy between the different determinations of [^f j^^^*^ can indeed be ascribed to higher 
order chiral corrections. 
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[^+]bare 




2 


5.5(2.7) - 4.9(2.4) 


0.56(14) - 0.59(14) 


2.5(1.3) 


3 


3.45(54) - 3.10(48) 


0.59(7) - 0.62(7) 


1.73(35) 


4 


2.75(31) - 2.46(28) 


0.62(5) - 0.65(5) 


1.47(23) 


5 


2.62(13) - 2.36(12) 


0.66(4) - 0.69(4) 


1.47(12) 


w.a. (|!^| > 2) 


2.68(11) - 2.41(11) 


0.63(3) - 0.66(3) 


1.49(10) 




1.1 


0.45 


0.25 



Table 10: Values of the bare couplings obtained from a NLO fit in the interval \zq\ < a for 
lattice Ai and xq = 7a, yo = 9a. The ranges indicated correspond to considering [FVf = 
1.14-1.40. 

4. Conclusions 

The purpose of this paper has been to estimate the weak low-energy couplings gf, defined 
in eq. (2.5), in the SU(4) chiral limit. Our method has been to measure the topological zero- 
mode contributions to three-point correlation functions of two pseudoscalar densities and a 
weak operator in sectors of non-trivial topology. The results of the measurements have been 
matched to NLO predictions of e-regime chiral perturbation theory. 

We have considered several fitting strategies for estimating the couplings gf and g^ , in 
an attempt to quantify the uncertainty induced by unknown higher order chiral corrections 
(NNLO), which are expected to be significant at the volume we have considered. While we 
observe small variations in the determination of g^ , as well as in the product gfgi, between 
the different methods, the value of g^ seems to be significantly affected by higher orders. 

Taking the bare couplings [fi'j^]'^^'^'' and [gigi]^^'^'^ cited in Table 8 for the lattice Ai, and the 
same renormalization factors and Wilson coefficients that were used in ref. [17], and inserting 
everything into eq. (3.10), we obtain 

g+ ~ 0.46(5) , g+g^ ~ 1.2(2) . (4.1) 

The errors shown involve statistical uncertainties as well as the uncertainty from the de- 
termination of the (quenched) pion decay constant F. We have checked that discretization 
effects in these numbers are small. However, systematic errors related to higher order chiral 
corrections as well as the quenched approximation have not been quantified. In any case, 
both quantities are in good agreement with those of ref. [17], where they were extracted from 
observables with very different chiral corrections, so this is a strong indication that higher 
order chiral corrections could be under control. 

The situation with g^ is less clear. The different fitting strategies we have explored give 
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values that differ by up to 30%, which is also the naive expectation for the magnitude of 
higher order chiral corrections. We have argued that an indirect extraction of from the 
combinations in cq. (4.1) is the method that should be least sensitive to the uncertainty 
induced by higher orders. With this approach, we obtain from the results of Table 9 

9T - 2.8(4) , (4.2) 

which is also in good agreement with the result of rcf. [17]. We should stress however that 
we have seen evidence that higher order corrections could be significant for this quantity; the 
corresponding systematic error cannot be quantified precisely, and has not been included in 
eq. (4.2). Probably an error of 30% in eq. (4.2) would be a reasonable estimate. 

The conclusions concerning the AI = 1/2 rule are the same as in ref. [17]: there is a 
significant AI =1/2 enhancement already in the SU(4) limit, which cannot be explained by 
penguin dominance. However, this enhancement is not as large as the experimental one. 

The method of the present study and that in refs. [16, 17] can be compared on two accounts. 
First of all, there is the issue of how well chiral perturbation theory converges with a given 
box size and geometry. In a symmetric box of size {FVf = 1.1, for instance, the magnitude 
of next-to-leading order corrections in the method of refs. [16, 17] is 15% (cf. Fig. 2 of 
ref. [16]), while in the present observables it is ~ 60% (cf. Table 1). In the present case 
the magnitude of the corrections can be reduced very significantly by taking the sources 
further away from the weak operator, while in the method of refs. [16, 17] this has no effect. 
Nevertheless, it could be concluded that from the point of view of chiral perturbation theory, 
the method of refs. [16, 17] appears to be preferable. 

The second comparison concerns the numerical cost of the measurements carried out. On 
this account, on the contrary, the present method appears to be preferable: a good statistical 
signal could be achieved with significantly less computational effort than in ref. [17]. Indeed, 
the number of quark propagators required per quark mass to construct the observables is a 
factor of 5 smaller in the present work. Note also that this factor scales with Ni^^, which is 
expected to scale with the volume^. 

These two competing aspects probably mean that, moving towards SU(3) symmetry and 
unquenched simulations, it would be wise to continue to probe the weak low-energy constants 
with (at least) two independent methods. In particular, the fact that chiral corrections are 
very different in the two cases, offers a good way of quantifying the systematics associated 
with the chiral fits. On the other hand, particularly in the SU(3) case where the penguin 
contractions need to be evaluated, which entails a significant numerical cost, it may be that 
the method introduced in the present work becomes preferable. 

®A new method to solve the F'^-problem of low- mode preconditioning has been presented in ref. [34], and 
in principle could also be applied to low-mode averaging. 
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Appendix A. Graph-by-graph results for the two-point correlator 

For completeness, we present in this appendix graph-by-graph results for the two-point cor- 
relation function defined in eq. (2.22). 

As usual, the Goldstone field is factorised into non-zero and zero-mode parts: 



Uix) = U^{x) Uo , U^{x) = exp 



F 

The propagator of the non-zero modes ^, which are perturbative, is of the form 



(A.l) 



iCijix) Ckiiy)) = 7; Sii5jkG{x -y)- 5ij6kiE{x - y) . (A.2) 



2 L 

Here G{x) is the massless non-zero mode propagator 

1 ^ /, j4)\e'^P'^ , , 27rnn 27rn. 



while E{x) is the "trace part" whose form is affected by quenching; in the unquenched case 
it reads E{x) = G{x)/Ni, while in the quenched case. 



^ 7TwG{x) + ^F{x) , (A.4) 

where 



2A^c 2A^c 



Since the zero-mode field Uq is an x-independent constant, and we are only interested in 
contributions to the correlation function {T"^{x)L\f{z))u that remain non-zero after taking the 
derivative dx^, it is clear that the only graphs that can contribute are those where the two 
operators are connected by a non-zero mode propagator. Denoting the operator V"" by an 
open square; the operator >Cq by an open half circle; the propagator in eq. (A.2) by a solid 
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line; the mass term in the chiral Lagrangian by a closed circle; the mass term originating 
from the Haar measure by a cross; /i = rriY^V; and choosing to list the results before taking 
the spatial average, time derivative, and zero-mass limit in eq. (2.22), we are led to: 



O 

□ Q 



□ — X— a 



n. 



a 



□ Q 



y9oG(z-x)(Tr {UoT^ + T''U^)T'' 



(A.6) 



'4P2 



[iVfG(O) - E{0)]doG{z - x) X 
Tr [(C/oT" + T''U^Q)T^]Tr [Uq + U^] 

Tr [{UoT" + T''U^Q)T^]j(Tr [Uq + U^q 
iNiT, 



3F2y 

imT?' 
x(^Tr 



00^(2 - s)G{s - x)(Tr [(C/qT" + T''Ul)T^^ 



(A.7) 
(A.8) 



doG{z — s)G{s — x) X 



+ 



2F2 
x<^Tr 



([/o^'^ + ^"c/o^){c/o + [/o^^^} 

5oG(z — s)£'(s — x) X 



6^2 L 



^JoT" + T"[/o^j Tr \t\Uo + C/(J) 
a2G(0) / doG{z - s)G{s - x) 



(A.9) 
(A.IO) 



+G(0) / aoa^G(z - s)a^G(s - x) (Xr {UoT'' + T''Ul)T^ 



0, 



|P^G(0)5oG(z - x)^Tr ^(C/qT'^ + r«C/o^)r^ 



(A.U) 



(A.12) 



~3F2 
x/Tr 



NiG{<d) - -E{0) doG{z - x) X 



(A.13) 



In the quenched case, there are two additional contributions, discussed around eqs. (2.27), 
(2.28): 



iC-term 
□ 



2iiyKNc 
mlFV 



doG{z - x) X (Tr {UqT'' - T''ul)T^ ) , (A.U) 
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A"- weight, conn. 



doG{z — x) X 



Tt 
-(TV 



(C/qT" + T''U^q)T^] Tr {Uo - U} 



{UoT" + r"C/T)r^ ) ( Tr {Uo - U^o 



(A.15) 



Here we made use of the fact that the integral over the zero-mode of the singlet field, ^o, is 
Gaussian, and $o can be approximated by the corresponding saddle point value. 



-2i 



mlFV 



(A.16) 



In fact this value was already inserted in order to arrive at the iT-term of eq. (2.28). 

As far as the zero-mode integrals appearing in eqs. (A.6)-(A.13) are concerned, the general 
trick to use is that, because of the invariance of the integration measure. 



{Aij Bki) = Cl 5ij5ki + C2 5ii5ji 



(A.17) 



where A,B ^ {Uq, Uq}. Carrying out contractions and solving the quadratic system yields 



Cl 



C2 



—^±—^i^Nf{TT[AB]) - {Tt[A] Tv[B])] . 



(A.18) 
(A.19) 



This leads to integrals whose values are listed in appendix B of ref. [21]. For the zero- mass 
limit in eq. (2.22) we only need the poles in l//i", which are also listed in ref. [21]. 
Applying this recipe in practice, we obtain 



lim {mVTj)^ 

m— >0 



lim (myS)( Tr {UoT"- + T^-U^T" 

m-+0 



Tr [{UoT^ + T'^C/(J)r*]TV [Uo + U^q] 



2|z^|Tr[rT^] 



Tr [{UqT'' + T^Wq)T^] ) (Tr [Uo + U^q] 



-4|zy|Tr [T'^r 



(A.20) 



(A.21) 



lim (myS)"( Tr {UqT'' + T'^;7o^){;7o + U^q^T"} 



4|z^|(2|i/| - iVf )Tr [T'^T^] , (A.22 



lim (mVT,) ( Tr 



UoT" + T^'UI 



Tr 



T'iUo + U^o] 



arpbi 



-4|z^|Tr [T"T' 



The additional integrals needed in the quenched case (eqs. (A. 14), (A.15)) read 



lim [mVT,)^ 

m,~>0 



Tr [([/qT"' + T''uI)T^]Tt [Uo - Ul 
Tr [{UoT" + r";7o^)r^])(Tr [Uo - C/j 



-2z^Tr [TT^] 



4i/ Tr [T^T 



(A.23) 
(A.24) 

(A.25) 
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Inserting the last two into eqs. (A. 14), (A. 15), we see immediately that the two quenched 
terms cancel against each other. 

Let us finally consider the spacetime dependence. After taking the spatial average and 
time derivative in eq. (2.22), omitting contact terms, and denoting Tx = {xq — zq)/T, we get 



dxo / dQG{z - x) 

J X 

doG{z - s)G[s - x) 
dxo / / doduG{z - s)dyG{s - x) 

Jx. J s 



1 

~T ' 
Thi{Tx) 
1 

T ' 



(A.26) 
(A.27) 
(A.28) 



where the function /ii(t) is given in eq. (2.26). In dimensional regularization, the object 
(?^G(0) appearing in eq. (A. 10) evaluates to 3^G(0) = 1/V . On the other hand, the constants 
G(0),i?(0), appearing in several contributions, cancel completely in the final result. 

Summing now all the results together, but making no assumptions about the form of E{x), 
and expressing the result as in eq. (2.22), we obtain 



B„{xo - zo) = — 



2|i/|T2 ^ ^ 2T ^ 



XO 



doG{z - s)E{s - x) 



(A.29) 



Appendix B. Graph-by-graph results for the three-point correlator 

For completeness, we present in this appendix graph-by-graph results for the three-point 
correlation function defined in eq. (2.30). To be precise, we list results for the operator Orsuv 
from eq. (2.32), with r, s, u, v assumed to be all different; results for the operator Oi are then 
obtained by symmetrizing according to eq. (2.31). 

Since the zero-mode field Uq is an x-independent constant, and we are only interested in 
contributions to the correlation function —{V°'{x)Orsuv{z)V^{y))u that remain non-zero after 
taking the derivatives dxgdyQ, it is clear that the only graphs that can contribute are those 
where the pseudoscalar densities are connected to each other or to the weak operator by 
non-zero mode propagators. Furthermore, in the SU(4) limit we can assume all the indices 
r,s,u,v to be different, which allows us to omit structures like 5ur, ^us, ^vr, ^vs', this means 
that the weak operator needs to be connected to at least one of the pseudoscalar densities. 
Denoting the operator "P" by an open square; the operator Orsuv by an open circle; the 
propagator in eq. (A. 2) by a solid line; and choosing to list the results before taking the 
spatial averages, time derivatives, and zero-mass limit in eq. (2.30), we are led to: 

□ OB = —df,G{x - z)df,G{y - z) x 



(B.l) 
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□ — H—O—O 



□ • O □ 



□ O o □ 



8F2 



[NfGiO) - E{0)]d^Gix - z)d^Giy - z) x 



+ {a ^ b) 



(B.2) 



d^Gix - z) jd^G{y - s)G{s - z) + {x ^ y) 



X <^ d^G{x - z) / d^G{y - s)G{s - z) x 



(B.3) 



8F2 



+df,G{y - z) / d^G{x - s)G{s - z) x 

J s 

x(^{UoT'' + T^U^o)ur{Uo + ul C/oT" + T'^f/jJ} 
~2d^G{x -z) [ d^E{y - s)G{s - z) x 

J s 

x{(UoT'' + T^U^uriUo + U^^sTt {UoT' + T'Uo) 
~2d^G{y -z) I d^E{x - s)G{s - z) x 

J s 

x{(UoT' + T'uluriUo + U^^sTt {UqT^ + T'^UI) 
+ (r <-> s, u <-> , 



(B.4) 



12F2 



G(0)a^G(x -z) d^d,G{y - s)d,G{s - z) 
+dlG{{))d^G{x -z) d^G{y - s)G{s - z) + (x ^ y) 
{UqT^ + T'^U^uriUoT'' + T'U^^s) + (a ^ 6) 

' G{x - s)G{y - s)d^d,Gis - z)d^d,G{s - z) 



(B.5) 



12F2 



+d,G{x - s)G{y - s)d^G{s - z)d^d,G{s - z) 
+G{x - s)d^G{y - s)d^,G{s - z)d^,dyG{s - z) 
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a 



-o— □ 





+d^G{x - s)d^G{y - s)df,G{s - z)d^G{s - z) 



{UoT'' + T''ul)UUoT' + T^f/o^).. ) + (a ^ 6) 



3F2 



-E{0) - NfG{0) 



d^G{x - z)d^G{y - z) x 



S2 r r i2 



Y^^G{x -y)^df,G{y - z)^ +{x^y)jx 
{UoT'' + T^U^o)us{UoT' + T'ulUr) + {a^b) 



Y? 



dlG{^)G{x - z)G{y - z) 



■G{^)d^G{x - z)d^G{y -z)x 



In the quenched case, there are two additional contributions, 



X-term 



-^G{x - y)d^G{x - z)df,G{y - z) x 



{UoT^ - T^ul)UUoT^ - T'U^o)vr ) + (a ^ 6) 



+ ^-^(^ ~ - z)d^G{y - z) x 

0, 
0, 



12^2 

+3G{0)d^G{x - z)d^G{y - z) 



(UoT'^ + T^U^o)usiUoT' + T''U^o)„. ) + (a ^ 6) 



(C/oT'^ + T''U^o)ur{UoT'' + T'U^oU) + {a ^ b) 
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■bTTU 



weight, conn. 



+ {a ^ b) 

umY?KNc 
2mlF 



(B.13) 



df,G{x - z)d^G{y - z) x 



(C/oT" + T^U^o)ur{UoT' + T'U^),,)(tt {Uq - u} 



+ (a ^ b) 



(B.14) 



Employing the same trick as in eqs. (A.17)-(A.19), the zero-mode integrals become 



limJmyS)2<| ( (C/oT" + r"[/t)„,(C/oT^ + T^ul),s ) + {a ^ b) 



2rp{arpb} _ U,\rr{a.rpb} 



limJmyS)2<j ( (C/oT" - r«[/t)„,(c/or^ - r^c/o^),, ) + (a ^ 6) 



(B.15) 



(B.16) 



hm (myS)3<j ( (f/oT" + T'^ulUiUoT^' + T^C/j),, )(Tr (C/q + Ul 



iUoT'' + r"C/ot)„,(C/or^ + T'U^oUTt {Uo + f/o^)) + (a ^ 6) 



16 



(B.17) 



limJmyS)3<| ( ([/oT" + r'^C/o^)„,{C/o + U^T^ + T^f/o^},, ) + (r ^ s, n ^ t;] 



(1^1 _ N,v^ + 2i.Vl)rM + - ^^')TtT^l 

2rp{arpb} , UAT^arpb} 

-^ur ^vs ' \ us vr 



(B.18) 



lim {mVT.f{ ( (?7or" + r'^?7j)„,(C/o + Ul^Ti {UqT' + T^f/J) ) + (r ^ s, u ^ f ) 

m — ^0 



(B.19) 



where W6 again omittGci all terms proportional to Sus^Sui'^i^vs^i^vr' 

The additional zero-mode integrals needed in the quenched case read 



m— >0 ' ^ 
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-8iy 



I |T-.{aj.&} _ rp{arpb} 
\ \ ur vs us vr 



(B.20) 



liin (myS)3|^(C/oT'^ + T''ul)uriUoT'' + T^U^UTt {Uq - ul) 

-((C/oT" + T'^ulUiUoT^ + r^C/(J)..)(Tr (C/q - ul)) + (a ^ 6)| 



= 16z/ 



I I tir t)s us 



(B.21) 



Inserting the last two into eqs. (B.13), (B.14), we see immediately that terms proportional 

to K cancel against each other. 

Let us finally consider the spacctime dependence. After taking the spatial averages and 
time derivatives in eq. (2.30), omitting contact terms, and denoting Tx = {xq — zo)/T, Ty = 
(yo - zo)/T,we get 



dxo / G{x - z) = h[{Tx) , 



yo 



d^G{y - s)G{s -z) = -S^oThi{Ty) , 



y -Js 



[s - z 



(B.22) 
(B.23) 
(B.24) 

(B.25) 



dyo [ [ d^dMy-s)dM 

Jy J s 

Oxodyo / G{x - y)\d^G{y - z) 

= y{-fi(.ry) + 2h[{Tx - Ty) [h\{Ty) + 2 IpfCMCpiTy)] } , (B.26) 

dxodyo / G{x - y)d,^G{x - z)dnG{y - z) 

= y{-h'l{Tx)h'l{Ty) + h'l{Tx - Ty) h'l{Tx) - h'l{Ty) + hl{Tx - Ty)^ (B.27) 

^ (B.28) 



= -yH{Tx,Ty) , 



8 8 



G{x - s)G{y - s)8^d^G{s - z)df,duG{s - z) 

+8^G{x - s)G{y - s)d^G{s - z)8^d^G{s - z) 
+G{x - s)d^G{y - s)d^G{s - z)d^d^G{s - z) 

+d^G{x - s)d^G{y - s)df,G{s - z)d^,G{s - z) 

= + ^{hli^x)h'l{Ty) + 7h'l{Tx - Ty) hl{Tx) - h[{Ty) 
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-3[/ii(r^ - Ty) + /ii(r^) + /ii(r^)] - 4[/i(r^) + /i(rj,)] } 



rx-Ty)\p\ 



CAryKiry) - Cp(r,)C;(r,)] } . (B.29) 



The functions hi, fi, Cp appearing here have been defined in eqs. (2.26), (2.36), (2.37), 
respectively, and in eq. (B.28) we identified the function H defined in eq. (2.35). We also 
need to know that in dimensional regularization, d^G{0) = 1/V and 

Pi 



G(0) 



V 



(B.30) 



where /?i is a "shape coefficient" [28, 29]. It is furthermore useful to note the identity 



) h[{T,) - h[{Ty) = hi{T, - Ty) + hi{T,) + hiiTy) . (B.31) 



Summing all the results together, but making no assumptions about the form of E{x), we 
obtain 



lim imVf T^d,,dy, [ [ {-V\x)Orsuv{z)V\y)) 



2rp{arpb} _ |,,|7^{a7'&} 
ur vs \ \ us vr 

2T 



X U + 



2|z/|r2 



+ 



doG{z - s)E{s -x) + dya 



doG{z - s)E{s - y) 



y -Js 



F2 



dxodyo I E{x - y)df,G{x - z)d^G{y - z) \ 



+ 



IjAl^lajib} I 2rp{arpb} 
\ \ ur vs ' us vr 

rp2 r 



2G{0) 



F^V 



hi^x) + h{Ty) - hl{Tx) - hl{Ty) + H{Tx, Ty] 



(B.32) 



Taking finally the combinations in eq. (2.31); writing the result in the form of eq. (2.33); 
identifying expressions of the form in eq. (A. 29) from the result; and inserting eq. (B.30) as 
well as the definition p = T/L, we obtain 



Aiy{xo - Zo,yo - Zq) ± Auixo - ZQ, yo - Zq) 

= (l =F |'S,,(xo - zo)B^{yo - zq) 



+^5 d 

^ p2 ^0 2/0 



/ E{x - y)d^G{x - z)d^G{y - z) 



FW 



2/5ip-2 + /i(r^) + /i(r^) - /^^(r^) - hi{Ty) + H{T:,,Ty) ) . (B.33) 
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